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MODELING AND SIMULATION OF 
THERMO-FLUID-ELECTROCHEMICAL ION FLOW IN 
BIOLOGICAL CHANNELS 
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Abstract. In this article we address the study of ion charge transport in the 
biological channels separating the intra and extracellular regions of a cell. The 
focus of the investigation is devoted to including thermal driving forces in the well- 
known velocity-extended Poisson-Nernst-Planck (vPNP) electrodiffusion model. 
Two extensions of the vPNP system are proposed: the velocity-extended Thermo- 
Hydrodynamic model (vTHD) and the velocity-extended Electro-Thermal model 
(vET). Both formulations are based on the principles of conservation of mass, 
momentum and energy, and collapse into the vPNP model under thermodynamical 
equilibrium conditions. Upon introducing a suitable one-dimensional geometrical 
representation of the channel, we discuss appropriate boundary conditions that 
depend only on effectively accessible measurable quantities. Then, we describe 
the novel models, the solution map used to iteratively solve them, and the mixed- 
hybrid flux-conservative stabilized finite element scheme used to discretize the 
linearized equations. Finally, we successfully apply our computational algorithms 
to the simulation of two different realistic biological channels: 1) the Gramicidin-A 
channel considered in US]; and 2) the bipolar nanofluidic diode considered in [30] , 


1. Introduction and Motivation 

Ion channels are ubiquitous in the cells of the human body. They allow communi¬ 
cation between the intra and extra-cellular sites and are responsible for the signaling 
pathways that regulate the activity of each part of the body system. A complete 
presentation of ionic channels is far beyond the scope of this article. For a treatment 
of the general properties of ionic channels and of their characterization in cellular 
biology and neuroscience we refer to the books Hg and Pi In broad terms, ion 
channel are “biological gates” that can be activated (opened) or deactivated (closed) 
by the application of electrical, chemical, mechanical or thermal stimuli. 

In the present article we mainly focus on the interaction among electrical, chemical 
and thermal driving forces, starting from the evidence that the sensory system in 
mammals is capable of discriminating thermal stimuli ranging from noxious cold (< 
8° C) to noxious heat (> 52° C). Of particular relevance are the biological temperature 
sensors denoted thermo-transient receptor potential (thermo-TRP) cation channels. 
Thermo-TRP channels constitute a “superfamily” of receptors in sensory neurons 
that can be activated by noxious stimuli, resulting in the transmission to the spinal 
cord and brain of the information of local pain perception [Hj . The fact that the 
thermal thresholds of many thermo-TRP channels can be modulated by extracellular 
mediators has recently led to the development of antagonists for noxious channel 
blockage in therapeutic uses as novel analgesics [32]. This is the case, for example, 
of heat-sensitive TRP channels exposed to capsaicin, a natural ingredient of spicy 
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foods such as red hot chilli peppers and of cold-sensitive TRP channels exposed to 
menthol. In [7] and m, it is shown how the chemical agonists capsaicin and menthol 
function for both types of TRP channels as gating modifiers, shifting activation 
curves towards physiological membrane potentials. 

Despite a significant amount of experimental data, the mechanisms underlying the 
marked temperature sensitivity of channel gating are still largely unknown. For this 
reason, the use of mathematical tools is increasingly becoming popular to support 
biophysical conjectures and suggest novel theories for data interpretation. 

In this context, the most widely adopted approach to ion channel modeling and 
simulation is represented by compartmental models described by a system of ordinary 
differential equations (ODEs) based on the solution of Kirchhoff current law (KCL) 
written at the cellular membrane level m The KCL equation is supplemented by 
phenomenological expressions characterizing the input-output functional response 
of each protein channel to changes in ion concentrations and/or electric potential 
inside and outside the cell [113] • Temperature in these models is usually assumed to 
be a given parameter. 

More sophisticated approaches involve the solution of a system of partial differential 
equations (PDEs) expressing balance of mass of each single ion species flowing across 
the channel and the Gauss law for the electric field. The system is supplemented by 
a transport relation, known as the Nernst-Planck (NP) equation, that describes ion 
motion under an electrochemical gradient mmm Also in this kind of modeling, 
well-known as the Poisson-Nernst-Planck (PNP) system and as the Drift-Diffusion 
(DD) system [25, T8], temperature is a given parameter. 

A significant step forward to account for temperature as a dependent variable 
was taken in DU- In this reference a hydrodynamic (HD) formulation including 
convective and thermal energy in the electro-chemical motion of a single cation 
is proposed and numerically investigated. A remarkable feature of m is that 
model and analysis are inspired and guided by the analogy between a biological 
ion channel and the channel of a semiconductor device in which electrons and 
holes, instead of charged ion particles, flow to transport electrical current between 
device terminals. This similarity between biology and solid-state electronics has 
been thoroughly addressed in the overview paper [12], in the numerical simulations 
of [15], 17) [2S] and in the MSc thesis [123] - In this latter reference, the mathematical 
view of m is generalized by the introduction of a hierarchical modeling perspective 
to represent ion transport in a biological channel in which the interstitial electrolyte 
fluid is assimilated to the semiconductor device medium where two monovalent 
species (anion and cation) are flowing under the effect of electric, chemical and 
thermal forces. The hierarchy proposed in |23] is based on the ideas discussed in 
the semiconductor modeling reference books [3B], [25], [IB] and [21], and includes 
four members: the basic DD formulation, the electro-thermal (ET) formulation, the 
hydrodynamic (HD) formulation and the thermo-hydrodynamic (THD) formulation. 
For each member of the hierarchy, a set of conservation laws for mass, momentum 
and energy is written. In the case of DD, ET and HD models, the energy exchange 
between particles and medium are neglected while in the case of the THD model 
a supplementary conservation law is added to account for the dynamical thermo¬ 
electrochemical balance among the three interacting subsystems. The hierarchy 
is extensively investigated in a series of numerical computations performed in a 
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simplified one-dimensional channel geometry. Externally applied data on which a 
sensitivity analysis is carried out are the values of bulk ion concentrations in the 
intra and extra-cellular sites and the applied potential drop across the channel. 
Simulations indicate that: i) channel heating is, in general, relatively small compared 
to ion heating; and ii) for certain ranges of model parameters, the high-order effects 
introduced by the THD picture can significantly affect the input-output transfer 
characteristics of the “biological transistor”. 


Based on the experience and results described above, in the present article we 
propose the following mathematical structure for the modeling of ion charge transport 
in biological channels: 


(a) : in Section [2] we review the classic DD (PNP) model for ion electrodiffusion; 

(b) : in Sections] we use the channel conformation model proposed in [2] to 
supply the DD (PNP) formulation with a set of “reduced order” boundary 
conditions that allow one to account for the experimentally accessible values 
of applied bias, ion concentrations and bathing temperatures; 

(c) : in Section[4]we review two members of the hierarchy of models for ion charge 
transport introduced in [23], the Thermo-Hydro dynamic and Electro-Thermal 
systems, and include electroosmotic effects by adding a linear advective term 
into the momentum balance equations proportional to the given electrolyte 
fluid velocity. The resulting modified hierarchy is a family of (approximated) 
velocity-extended models, as the self-consistent study of electrolyte fluid 
dynamics through Navier-Stokes equations is neglected unlike in (HT, THJ [201, 

E5j; 

(d) : in Section [5] we illustrate the solution map used to solve iteratively the 
models in (a), (b) and (c) in steady-state conditions; 

(e) : in Section [6] we study a linear model advective-diffusive and reactive 
boundary value problem (BVP) that represents each of the subproblems 
arising in the solution map, focusing on the continuous maximum principle; 

(f) : in Section [7] we address the numerical study of the model BVP in (e) based 
on the dual-mixed hybridized (DMH) finite element method proposed and 
analyzed in [31 13 12Sj for elliptic equations and extended to more general 
differential operators in [2]; 

(g) : in Section [8] we propose a stabilization of the DMH scheme to deal with 
the case where advect-ion or reaction dominate over diffusion in the model 
BVP; 

(h) : in Section [ 9 ] we conduct a series of numerical tests to verify the convergence 
and stability properties of the DMH method; 

(i) : in Section [lO] we apply external temperature gradients, besides the usual 
electrochemical gradients, to the biological transistor to replicate in vitro the 
biophysical situations occurring in vivo in thermo-TRP channels; 

(j) : in the concluding Section [Tl] we summarize the main contributions of the 
present research and we address a list of forthcoming activities for model 
improvement. 
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2. The Poisson-Nernst-Planck model for ion electrodiffusion 
The Poisson-Nernst-Planck system (PNP) for ion electrodiffusion reads |dlj : 


(la) 

Qc- 

-q l- + divfj (ci,(p) = 0 

i = 1, • 

■ ■ ,M 

(lb) 

f i ((k, V) = -DiWci + /h—CqE 

\%i \ 

i = 1,. 

■ ■ ,M 

(lc) 

Di = lllVth 
\Zi\ 

i = l,. 

■ ■ ,M 


(l d) divE = — 

(le) E = — V<p 


M 

(If) p = qJ2 z iCi + qP. 

i =1 


In the PNP system, (la) is the continuity equation describing mass conservation 
for each ion whose concentration is denoted by c* (m” 3 ), i = 1 M > 1 

being the number of ions flowing in the electrolyte fluid. Each ion flux density 


(m 


—2„—1i 


is dehned by the Nernst-Planck relation (lb) in which it is possible to 


recognize a chemical contribution and an electric contribution, in such a way that 
the model be regarded as an extension of Fick’s law of diffusion to the case where 
the diffusing particles are also moved by electrostatic forces with respect to the fluid. 
The quantity z % is the valence of the i-th ion, while /q and _D,; are the mobility and 
diffusivity of the chemical species, respectively, related by the Einstein relation (lc) 
where V t h = kBT sys /q (V) is the thermal potential, our having denoted by ks, T sys 
and q , the Boltzmann constant, the absolute temperature of the system and the 
elementary charge, respectively. We remark that in the PNP modeling approach, 
T sys is a parameter representing the constant temperature of the system, without 
distinction between the ionic species and the electrolyte fluid. This assumption is 
going to be relaxed in Section |4j The function p is the space charge density (Cm” 3 ) 
in the electrolyte and is given by the sum of two contributions, the mobile ion charge 
qJ2iL i z i c i and the fixed charge qP. The electric field E (Vrn” 1 ) due to space charge 
distribution p in the electrolyte is determined by the Poisson equation (Id) which 
represents Gauss’ law in differential form, e being the dielectric permittivity of the 
electrolyte fluid medium. For further development, it is useful to introduce the 
electrical current density jj (Am” 2 ), equal to the number of ion charges flowing 
through a given surface area per unit time and defined as 


(lg) j i : = qzdi = —qZiDj'Vci + qpi\zi\c-(E i = 1,..., M. 

Remark 1. The PNP system (|TJ) has the same format and structure as the Drift- 
Diffusion (DD) equations for semiconductors (see, e.g., jTHJy), but it is applied to 
a different medium (water instead of a semiconductor crystal lattice) and includes, 
in general, more charge carriers than just holes and electrons, as in the case of 
semiconductor device theory. 
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3. Geometry, biophysical assumptions and boundary conditions 

Following the presentation of [0], we illustrate in Fig.[l]the schematic representation 
of a cross-section of a biological channel, assuming rotational invariance around the 
channel axis (x axis). 



FIGURE 1. Schematic view of the biophysical problem. Five regions 
can be distinguished. From left to right: intracellular bathing solution, 
antichamber, channel region, antichamber, bathing solution in the 
extracellular side. At the endpoints, x = 0 and x — L, terminal 
contacts are marked in red color. 


Five regions can be distinguished, ordered from left to right: 

(1) Region 1, L < x < 0 _ : this is the bathing solution in the intracellular side. 

(2) Region 2, 0“ < x < 0: this is the channel antichamber (or access region) 
from the intracellular side. 

(3) Region 3, 0 < x < d: this is the channel region. 

(4) Region 4, d < x < d + : this is the channel antichamber from the extracellular 
side. 

(5) Region 5, d + < x < R: this is the bathing solution in the extracellular side. 

The above geometrical representation leads naturally to a multi-domain formulation 
for ion channel simulation. The biophysical treatment of such a problem is fully 
carried out in p] whilst its mathematical and numerical treatment is the object of 
the present article, where the PNP formulation (Jl]) is extended to include thermo- 
hydrodynamical phenomena to ion electrodiffusion in the channel. 

3.1. Assumptions. The endpoints of the domain, x = L and x = R, are located 
sufficiently far from the antichamber and channel regions, in such a way that 
appropriate equilibrium conditions can be applied. More importantly, the endpoints 
are assumed to be the physical place where the solution intra and extra-cellular 
electrochemical conditions are accessible to experimental measurements. Because of 
this, following M, we assume that at x = L and x = R: 













(2a) 

(2b) 

(2c) 

(2d) 
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(Al): the electric potential ip is a known given quantity, so that: 

<P(L) = <Pl 
<p(R) = ip R -, 

(A2): the ion concentrations c t are known given quantities, so that: 

Ci(L) = Ci : L, i — 1,..., M 
Ci(R) = c itRj i — 1,... ,M 

where M > 1 is the number of ions flowing in the cellular solution. The 
boundary values for the ion concentrations satisfy the clectroneutrality con¬ 


straint: 


(2e) 

(2f) 


M 

= 0 

i= 1 
M 

Z i°i,R = 0 

i= 1 


where Zi is the charge number associated with each ion species c t (zi > 0 for 
cations, Zi < 0 for anions and Zi = 0 for neutral species). 

We close the characterization of the bathing regions by assuming that: 

(A3): the ion flux densities vanish inside the baths 

(2g) Ji(x) — 0 L < x < CT, i — 1,..., M 

(2h) Jj(x) = 0 d + < x < R, i — 1,..., M. 

Continuing to move from the periphery of the domain towards the channel region, 
we encounter the antichamber openings, at x = 0~ and x = d + , respectively. Located 
on the membrane lipid bilayer, a space-dependent fixed charge density P = P(x) is 
distributed. We make the following assumptions within the antichamber regions: 
(A4): electroneutrality holds at the channel mouth entrances: 

(2i) qP( 0“) + q^ZjCji 0“) = 0 

(2j) qP(d + ) + qJ2zjCj(d + ) = 0; 

(A5): the electric potential is constant in the antichamber regions: 

(2k) ip(x) = <p( 0 _ ) 0“ < x < 0 

(21) <p(x) = ip(d + ) d~ < x < d + ; 

(A6): the ion concentrations are constant in the antichamber regions: 

(2m) Ci(x) = Cj(CT) 0“ < x < 0 i — 1,..., M 

(2n) Ci(x) = Ci(d + ) d < x < d + i — 1,..., M. 

3.2. Boundary Conditions at Channel Openings. In this section we use the 

geometrical multi-domain representation of the problem of Sect. [3] and the assump¬ 
tions made in Sect. |3.1| to derive the boundary conditions to be supplied to the PNP 
Model at x = 0 and x = d (channel openings). Using Einstein’s relation (lc) in (lg) 
we can write the PNP current density as 

d^V 


Ji Q \%i ■ 


dx 


(3a) 


i — 1,..., M 
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where 

(3b) 


fT '■=f + —V t h In 


Cj 

c c ref, 


i — 1,..., M 


is the electrochemical potential and c re f := max {cj Cj /?} is a reference concentra- 

tion. Inverting (3b) we obtain the well-known Maxwell-Boltzmann (MB) statistics 
for the ion densities 


(3c) 


Ci - c reJ exp I Zi 


PC 

Vi-V 


V th 


i = 1, 


, M. 


Using (A3) we get 

(3d) vTi x ) = const =: ip^ L L < x < CT 

(3e) vT( x ) = con st =: (pf R d + < x < R 


i — 1,..., M 
i = 1. M 


where tpf c L and <Pi C R , i = 1,, M, are constants yet to be determined. From (3c) we 
get 


(3f) Ci(x) — c ref exp ( Zi 


(Vl c (x) -ip(x)) ' 


L < x < R i — 1,..., M. 


Thus, using (3d) and (3e) and assumptions (A5) and (A6), we find: 

(Wl ~ V( x )] 


Ci(x ) = c re f exp ^ 

Ci(x) = Cref ex P Zi 


v th 

(Wr ~ v{ x )) 

v th 


L < x < 0 


i = 1, 


, M 


d < x < R i — 1,..., M. 


Using (A2) in the equations above at z = L and z = R, respectively, we obtain the 
boundary values for the electrochemical potential: 


(3g) 


(3b) 


Vl + 

v th , 

— In 


i = 1 ,...,M 

Zi 

V c ref J 


Vr + 

^ln( 

f Ci,R \ 

i = 1,..., M. 

Zi ' 

\ Cref J 



Remark 2. R is important to note that relations (3g) and (3h) allow one to express 
the electrochemical potentials at the contacts as a function of the sole accessible 
quantities. This is in contrast with m, page 4%, relations (3.8a)-(3.8b), where the 
values of the electrochemical potentials were imposed and the corresponding values of 
the electrical potential were computed using (3.14a)-(3.14b). 


From the previous discussion, we see that electroneutrality holds at the external 
contacts (x — L, x — R) and at the channel mouth entrances (x = CU, x = d + ), but 
not, in general, elsewhere. Therefore, it makes sense to introduce the voltage drops 
occurring in the bathing solution regions, the so-called built-in potentials: 

(4a) <pbi{ 0~) := <p(L) - <p( 0") 

(4b) (fbi(d + ) := v(d + ) ~ v(R)- 
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Using the charge neutrality conditions (2i) and (2j), the MB relations (3c) and the 
definitions Q , the built-in potentials can be determined by solving the two following 
nonlinear algebraic equations: 


(4c) 

(4d) 


P (° ) + 2^ z i c i,L exp I Zi — 

i =1 V 

M / 

P{d + ) + 5Z z i°iJi ex P 


v th ) 
Md + y 


= 0 


2=1 


V, 


= 0. 


th 


Remark 3. In the particular case of the KCl solution (z\ = +1, Z 2 = — 1) considered 
in [S], equations (4c) and (4d) can be explicitly solved, to yield: 


Vbi{ 0 ) = Vth In 


Vbi{d + ) = V th hr 


- P ( 0 -) + ^( P ( 0-)) 2 + 4 ( c 1 , l ) 2 ~ 

201^ 

' P(d + ) + \] (P(d + )) 2 + 4(c 1i r) 5 


2ci, 


R 


where we notice that ci t L = c 2i l and C\jt = c 2) _r because of the electroneutrality 
constraints (2e) and (2f). 

In order to complete the characterization of the boundary values for the dependent 
variables of the problem we need the value of the electric potential and of the ion 
concentrations at the two endpoints of the channel. The potential is determined 
using Q and (A5), which yield: 

(5a) <p(0) = <p(0 - ) = <p L - ^6i(0") 

(5b) ip{d) = ip(d + ) = (f R + <pbi(d + ). 


To determine the ion concentrations we use (|3f|), (|3g|) and (|4a|) at x = 0 to obtain 
(5c) Cj(0) = c re f exp ( z : 




Vt 


th 




= Ci,L exp Zi- 


t P\n (0 


v t 


i = 1, 


M 


th 


and (3f), (3h) and (4b) at x = d to obtain 
(5d) 


Ci(d) = c ref exp Zi 


(<Pi C (d) ~^(d)) \ 


V t 


th 




= c ijR exp -Zi 


<Pu(d + y 


V t 


i — 1,..., M. 


th 


Summarizing, the Dirichlet boundary conditions for the PNP system in the reduced 
channel domain = (0, d) are at x = 0: 


M 


(6a) 

(6b) 

(6c) 

(6d) 


P(0 ) + ^^c,:,Lexp (::, 6'(' : .' 1 ^ | - 0 

V(0) =tp L - <pn((T) 


V, 


th 


<P?{ 0) = VL + —Vth In ( — 

\ Cref , 


i = 1, 


, M 


q(0) = c hL exp z 


I^bi (o 


Vt 


i — 1 ,...,M 


th 
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and, at x 
(6e) 

(6f) 

(6g) 

(6h) 


= d: 

P(d + ) + J2 z i c i,R ex P [ z i Pl "^ = 0 

i= 1 V v th ) 

<p(d) = <p R + <fbi(d + ) 

<P?(d) =<Pr + 1 V th In ( i — 1, ... ,M 

\Cref J 

Ci(d) — c ijR exp l—Zi ——- 1 * = 1 

4. Extensions of the PNP model 


In this section we present two extensions of the PNP equation system illustrated 
in Sections [2] and [3j From now on we restrict our attention to the study of a 
channel in stationary conditions and to simplify the exposition, in close analogy with 
semiconductor device physics, we consider a binary mixture of monovalent anions 
and cations, i.e., such that their chemical valence is equal to ±1 as in the case of the 
KCl solution or the NaCl solution. We use the symbol p (positive) and n (negative) 
to refer to the concentrations of cations and anions, respectively. The symbol v is 
used to indicate either p or n. When the operators =F or ± are used, the upper sign 
always refers to n-type particles, the lower sign to p-type particles, and the value for 
v must be chosen accordingly. The quantities associated with the electrolyte medium 
are referred to by the symbol e, while x denotes henceforth the spatial coordinate. 

In our analysis we follow [23], the series of references mmm in the context of 
semiconductor device modeling, and the theory of mmmm- 

The first extension is denoted velocity-extended Thermo-Hydrodynamic model 
(vTHD) and the second extension is denoted velocity-extended Electro-Thermal 
model (vET). Both extensions include thermal driving forces and the translational 
contribution of electrolyte fluid flow, besides the usual electrochemical forces ac¬ 
counted for by the PNP equations. In the vTHD model the temperature of anions, 
cations and lattice are considered as distinct dependent variables, while in the vET 
model there is only one temperature, generally varying with the spatial coordinate x, 
to describe the global system. 

Finally, in the present article we assume the electrolyte fluid velocity v e to be a 
given function of x. This is a strong simplification which is going to be removed 
in a future publication. Moreover, to self-consistently account for electrodynamical 
effects, we adjoin to both vTHD and vET model the solution of the Poisson equation, 
written below in conservation form: 


(7a) 

(7b) 


OE q. 

& = h p -" +p) 

E = ~^. 
dx 


4.1. Thermodynamical Equilibrium. Let T p = T p (x), T n = T n (x) andT e = T e (x) 
denote the temperatures of cations, anions and electrolyte, respectively, with x G [0, d\. 
Let also T sys be a given value representing the constant temperature of the whole 
environment (ions + fluid). By definition, thermal equilibrium is the condition 
corresponding to applying no external electrical, mechanical, chemical and thermal 
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forces to the biophysical system. In such a case, system response is represented by 
the following values of the dependent variables: 


(8a) 

V v (x) 

= 0 

(8b) 

V e (x) 

= 0 

(8c) 

Tp( x ) 

= T n (x ) = T e (x) = T sys 


Eq. (|8aj) expresses the fact that the ion drift velocity is equal to zero, while Eq. (8b) 
expresses the fact that the fluid translational velocity is equal to zero. Eq. (8c) 
expresses the fact that the system settles in an isothermal condition. The three 
above conditions imply that: 


(8d) J v {x) = 0 

(8e) S v (x) = S e (x) = 0 


where J p and J n 


are anion and cation current flows, while S u and S e are ion and 


fluid energy fluxes. Using (3a) in (8d) tells us that at thermal equilibrium each ion 


electrochemical potential is constant in [0, d] 


(8f) 


<Pt c (x) = V 


ec 

v,eq 


where the constant value p e ^ eq is different for each ion and can be computed as 
detailed in Sect. 3.2 Model consistency requires that under thermal equilibrium any 
extension of the PNP system verifies all conditions ([8]). This fact is documented in 
the later sections where the vTHD and vET models are introduced. 


4.2. Velocity-extended Thermo-Hydrodynamic (vTHD) model. The sta¬ 
tionary velocity-extended thermo-hydrodynamic model in the one-dimensional setting 
consists of the following system of conservation laws (see £321 Chap. 2] and [131 141132]): 


(9a) 

(9b) 

(9c) 

(9d) 

(Be) 

(9f) 

(9g) 

(9b) 

(9i) 


lcV, 

q dx 


= 0 


d ( J, 


dx V v 


dv 


d 


Jv + — = ±qD u - - qi^vV— p =F 


k B T v 


dx 


dx 


=F qW e 


= EJ„ 


v 


[w v - wl q ) 


— ( W P - W p) + —iw n - W eq ) 


T, 


T n 


dSu 

dx 

dSe 

dx 


dT J 

S v = -Kv^T^{!Ou + k B T v ) 
dx q 

dT e 

S e Ke~r\ T N s V e W e 
dx 

1 , l2 3 

Wv = -jj n v\ v v\ + 2 k BT v 

1 , ,2 3 

w e = -m e \v e \ + -k B T e 


w eq = -k B T e . 
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In system (J9]) , v and T v represent ion density and temperature, N e and T e are 
electrolyte number density and temperature, while v u , J u , S v and S e represent the 
drift (or translational) velocity of each ion, the ion current density, the ion energy 
density flux and the electrolyte energy density flux, respectively. The drift velocities 
and the current densities are related through the phenomenological law 


(10a) 


Vu 


Ju 
=F — 
q v 


Eq. (9a) expresses conservation of ion density, eq. (9b) expresses conservation 
of ion momentum density, eqns. (9c) and (9d) express conservation of the energy 
densities for ions and electrolyte fluid, eqns. (9e) and (9f) are constitutive laws for 
ion and electrolyte energy density fluxes S v and S e , respectively. 

The quantities k u and K e are the thermoconductivity coefficients of ions and 
electrolyte, respectively, while r Pv and r Wv are the momentum and energy relaxation 
times, respectively. Phenomenological expressions for these coefficients are: 


(11a) 

(nb) 

(lie) 


T, 


Vu 


VPI'uf^OuJ'e 

qTu 


3 Houk B T u T e 

2 qVsat(Tu + T e ) 

3 Ho k^T e 

x- za 


+ 


'Vu 


2 ’ 


where fj,Q v is the low-field mobility and v sat is the saturation velocity. For a physical 
interpretation of the parameter v sat and its influence on system behavior, we refer 

to [in, 23] • 

The ion mobility 


(12a) 


l^u 


m v 


'Vu 


is related to the diffusion coefficient D u by the generalized Einstein relation 


(12b) 


k B T„ k B T u 

u l^v Tp v 

q rriu 


where m u is the ionic mass. 

Finally, w u and w e are the ion and electrolyte energies, in which we can distinguish 
a kinetic contribution and a thermal contribution. The quantities Wp q and w are 
the equilibrium energies of cations and anions. 


Proposition 1 (Consistency with respect to thermal equilibrium). The vTHD 
satisfies relations (J8|. 


Proof. Using (8a), (8b) and (8c) in (9b) and (12b) we get 

j , kaTgys dv dip 

J v = ±qpu -— qhuv^- = o 

q ox Ox 


J u = -qp v v-^ (=F V th In ( V/Vref) + <p) 


0. 


or, equivalently, 
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The above equation represents the particle flux equilibrium condition (8d) and can 


be written in the form (3a) upon defining the equilibrium electrochemical potential 


■■= V =F V th bi 


v 


V, 


ref, 


which corresponds to (3b) having set z n — — 1 and z p = +1. Having proved (8d) 


(and, equivalently, (8f)) and using (8c) in (9e) and (9f), we immediately obtain 


S n = S p = S e = 0 


that is the energy flux equilibrium condition (8e). We notice that at thermal 


equilibrium the ion and fluid energies coincide with their rest energy 

IPn ICp ^^B^sys- 


□ 


4.3. Velocity-extended Electro-Thermal (vET) model. This model differs 
quite significantly from the vTHD because of the following simplifying assump¬ 
tions: 

(HI): cations, anions and electrolyte are in local equilibrium, i.e. 

(13a) T n = T p = T e =: T{x). 

This situation occurs when the system is assumed to have had enough time 
to evolve and to reach steady state; 

(H2): the kinetic energy is small compared to the thermal energy, i.e. 


(13b) 


(13c) 


Thus we can write 


1 , ,2 3 

-m u \v u \ 2 <C - k B T. 


w v ~ -k B T 


and neglect the nonlinear inertial term in the momentum balance equa¬ 


tion (9b). 


Using (HI), and summing up the three resulting energy balance equations we 
obtain the following equations of the velocity-extended ET system: 


(14a) 

(14b) 

(14c) 

(14d) 

(14e) 


-cliv J u = 0 

q 


d 


J u = qn u uE ± ii v — (uk B T) =F quv e 


div S tot = E • (J n + J p ) - div 
Stot = ^ k B TN e v e 

R tot T R n T ^p‘ 


5 knT 


2 q 


( J p J n ) 


Remark 4. Looking at the (simplified) momentum conservation equation (14b), we 
see that it now furnishes an explicit expression for the current density as a function 
of concentration, electric potential and temperature. This allows one to adopt the 
basic ideas of the Scharfetter-Gummel discretization method [34], which is the most 
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widely used numerical scheme in contemporary semiconductor device simulation 
because of its remarkable stability and accuracy. 


Remark 5. Simple manipulations of (14b) show that it can be written in the 
following equivalent DD form 


(14f) 

where 

(14g) 


dv 

Jv = ±qD u — + qp v vE v 


E v = Et 


± 


d (knT N 




dx 


q 


The effective electric drift E u acting on each ion species can be interpreted as the 
linear superposition of the electric, fluid and thermal forces that drive ion motion in 
the electrolyte channel fluid. Relation (14f) can be interpreted as the formal limit 
of (9b) as the relaxation time t Pv —* 0. 

Proposition 2 (Consistency with respect to thermal equilibrium). The vET satisfies 
relations 


Proof. The proof is similar to that for the vTHD model and therefore it is omitted. □ 


5. Solution map for the vTHD system 


In this section we describe the solution map that is used to iteratively solve the 
vTHD and vET models. The description of the method is discussed in detail in the 
case of the vTHD system as a similar approach is used to treat also the vET model. 
The adopted solution map is an extension of the decoupled algorithm known as 
Gummel’s map, a functional tool widely employed in contemporary Drift-Diffusion 
simulation of semiconductor devices P31 EEl EE]. The Gummel solution map is 
basically a nonlinear block Gauss-Seidel iteration that allows one to subdivide the 
considered PDE model system into blocks of single equations to be solved separately 
and in sequence until convergence is achieved. 

To describe the method, we follow the idea proposed in [I] (cf. Eq. (25) of this 
latter reference), and we introduce the generalization of the Maxwell-Boltzmann 
statistics (|3c|) to the case of the vTHD system: 


(15a) 

(15b) 


n = 


Cref exp ( q- 


p = c ref exp q 


_<p-<Pn 

k B T n . 

Tp T 
k B T p , 


where, for notational brevity, we have omitted the superscript ec in the electrochemical 


potentials for anions and cations. The main difference between the definitions (15) 


and the corresponding (3c) valid in the case of the standard PNP model, is that the 
ion temperatures T n and T p are used instead of the system (constant) temperature 
T sys . Assuming local thermal equilibrium at the channel entrance and outlet, the 
temperature of the ions and of the fluid satisfy the following Dirichlet boundary 
conditions: 


(16a) T„"(0) = T“(0) = T e °(0) = T L 

(16b) T«(d) = T°(d) = T°(d) = T B . 
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We notice that and T ri are the externally accessible temperatures of the intra and 
extracellular baths, not necessarily being equal to the same value. This feature of 
the model is important when an external temperature gradient has to be enforced 
across the channel to model the heat biosensors described in Sect, [l] Conditions (16) 
imply that Vth = k B Ti,/q at x = 0 and V t h = k B T B /q at x = d. 

The Gunnnel algorithm for the iterative solution of the vTHD model starts with 


an initial guess 
and (16). Then, 
steps: 

( 1 ) 


Tp , T°, T e ° that satisfies the boundary conditions (J6]) 
:or k > 0 until convergence, the algorithm consists of the following 


(17c) 


solve with the damped Newton method (see [36], Chapt. 7) the nonlinear 
Poisson equation (NLP) resulting from substituting (15) into ([7]). This returns 
the updated potential tp k+1 ] 

solve the continuity and momentum balance equations ([9a])- (9b). This 
returns the updated concentrations p k+1 , n k+l ; 


(9e) . This returns the updated ion tempera- 


( 2 ) 

(3) 

(4) 

(5) update the electrochemical potentials by inverting (Jl5|) : 


solve the energy equations (9c) 
tures Tp +1 , T k+l ; 
solve the fluid energy equation (9d) 
temperature T e fc+1 ; 


(9f). This returns the updated fluid 


(17a) 

(17b) 

( 6 ) 


ip k n +l = + +l 


V k p +1 = V 


k B T k+1 ^ ( n 
q 


k+1 , kaT^/p 


+ 


q 



Cref 


check the convergence of the iteration by controlling whether the maximum 
absolute difference between two consecutive iterations k and k +1 is less than 
a prescribed tolerance 


max 11 n 
v eu 


k+l 


V 


< toll 


where toll is a given tolerance, U is the set of unknowns [<p, ip p , ip ni T p , T n , T e 
and II • 


is the L°° norm. If condition (17c) is satisfied, the algorithm stops. 


A deep analysis of the convergence properties and mathematical foundations of 
the Gummel decoupling algorithm for the stationary DD model of semiconductors 
can be found in [TEJ • Preconditioning strategies for improving the convergence rate 
of the map are investigated in [23j. 


6. The Advection-Diffusion-Reaction Model Problem 

As outlined in items (e)-(h) of the introduction, we are providing a detailed 
description of the algorithm, beginning in this section, and continuing through 
Section 9. 

All the linearized equations of the vTHD system ([9]) that need be solved during 
the iterative procedure described in Sect. [5] can be cast into the following general 
Boundary Value model Problem (BVP): 
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find u 

and J such that: 



dJ 


(18a) 

— b cu — g 
ox 

in 

(18b) 

T n du 

J = vu — U— 

ox 

in Ll 

(18c) 

u = u 

on dVt 


In (18), SI is the open interval (0 ,d), u and J denote the primal unknown and the 
associated flux, while c and g are the reaction and production terms such that (18a) 
can be regarded as a stationary conservation law for u. In view of the numerical 
treatment of (18) we assume D, v, c and g to be piecewise smooth functions in 
O, the diffusion coefficient D being a positive bounded function while c and g are 
nonnegative given functions. Dirichlet boundary conditions are expressed by the 
nonnegative function u : dfl —> M + . Let 


(19a) 


Cu \= 


dJ 

dx 


+ cu = 



dii 

-D — + vu 

OX J 


+ cu : S7 —> 


The following properties express important biophysical features of the solution of 
the linearized BVP (18). We refer the reader to [3D] for details and examples. 


Definition 1 (Inverse monotonicity). Let w G C 2 (Ll) D C°(f2). We say that C is 
inverse-monotone if the inequalities 


(19b) Cw{x) >0 Vx G 

(19c) w(x)>0 V x G dLl 


together imply that 

(19d) w(x) > 0 VxGfl 


Remark 6. Inverse monotonicity expresses the fact that the dependent variable of 
the problem, say, a concentration, a temperature or a mass density, cannot take 
negative values. 

Theorem 1 (Comparison principle). Suppose that there exists a function 0 G 
C 2 (Ll) n C°(Q) such that 

(19e) jCu(x) < C<f>(x) W x G Lt 

(19f) u(x) < f>(x) VxG dLl. 

Then, we have 

(19g) u(x) < cf>{x) VxGfi 


and we say that <p is a barrier function for u. 


Combining (19d) and (19g), we obtain the following result which is a very useful 
tool in the approximation process of the BVP (18). 


Theorem 2 (Continuous Maximum Principle). Suppose that C is inverse-monotone 
and that the comparison principle holds for a suitable barrier function 0. Then, 
setting := max0(x), we have 


(19h) 0 < u(x) < Mfj, Vx G 

and we say that u satisfies a continuous maximum principle (CMP). 
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7. Finite Element Approximation 


The content of the present methodological section is organized as follows. In 
Sect. 7.1 we introduce the dual-mixed weak formulation of (18) and in Sect. 7.2 its 
corresponding hybridized finite element approximation, denoted DMH method. Then, 
in Sections A3 and 7_A we illustrate how to reduce the computational complexity of 
the dual-mixed hybridized method by the use of static condensation, which leads to 
solving a linear algebraic system of the same structure as a standard nodal-based 
finite element scheme. 


7.1. Continuous formulation. In the dual mixed method, the variables u and J 


in (18) are treated as independent unknowns, each of which belongs to a different 
functional space, namely u G V := L 2 (12) and J G Q H( div, 12) = 7/ 1 (fl). We pro¬ 
ceed by formally multiplying equation (18a) by a test function w G V and equation 


(18b) by a test function q G Q, and integrate by parts over 12 to obtain the problem: 


find (u, J) EV x Q such that, for all (w, q) G V x Q: 


(20a) 

(20b) 


f D 1 (J — vu) qdx — [ u-^— dx = — [ 
Jn Jn ox J do 


[ d /wdx 
Jn ox 


/ cuw dx 
Jn 


/ gw dx. 
Jn 


(q ■ n )u da 


Existence and uniqueness of the solution pair (it, J ) of problem (20) can be proved 
under suitable coercivity assumptions, see mi. 


7.2. Discrete formulation. We denote by 12 the open interval (0, L) and by <912 = 
{0,L} its boundary on which an outward unit normal n is defined, in such a way 
that n(0) = — 1 and n(L) = +1 We introduce a partition (triangulation) Th of 12 
into a number N d > 2 of intervals K l = [xi,x i+ i),i = 1 ,..., N e i, with x\ — 0 and 
xn b1 +i = L. On each element K t we introduce a local unit outward normal vector 
n| dKii w e denote by hj := x l+ \ — Xi the length of the interval and set h := max^ hi. 
We also set N := N e i + 1 and denote with A4 := the set of vertices of Th- 

We associate with Th the following function spaces: 

(21a) Vh := {wh € L 2 (12) s.t. Wh\x £ lPo(2 t)VA e Th} 

(21b) Q h := {q h e L 2 { 12) s.t. q h \ K G Pi (AT) VAT G %} 

(21c) A hiP := {X h : J\f h ->• R s.t. A h (a) = p a , A , h (b) = p b } 


where Vh. is the vector space of piecewise constant polynomials defined over 12, Q h is 
the vector space of piecewise linear polynomials discontinuous over 12, and A h, P is 
the space of functions defined only at the vertices of Th with given boundary values. 
Notice that dirn(Qh) = 2 N e i because its functions are not continuous at interelement 
interfaces. The DMH formulation of problem (20) is: 
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find {J h ,u h , \ h ) e Q h X 14 x such that, for all (q h ,w h ,^ h ) £ Q h xV h x A hfi : 

dq h 


(22a) 

(22b) 


E 

KeT h L' 


E 

AeT h 


( Db 1 (J h -vu h )q h clx - [ u h -^-dx+ f (q h ■ n\ dK )X h dx 
J K J K Cs 'X J dK 

f dJ h 


= 0 


Ik dx 


Wh dx 


IK 


cuhWh dx 


= - E / 9 w h dx 


K&Th ' 


I< 


(22c) / {Jh ■ rv\a K )ihdx = 0 

J dK 


Ker h 

where Dh : Th —> 1R + is a modified diffusion coefficient containing a stabilization 
term to deal with advective-dominated problems. 


Remark 7. The proof of existence and uniqueness of the solution ( Jh,Uh , A h) of (22) 


is not an easy task because of the advective term. In what follows, we assume that 
such a solution always exists and we refer to I5U2H] for the general theoretical tools 
required in the analysis of (22) and to |2] for an example of such analysis in the case 
of mixed methods interpreted as nonconforming discretizations. 


Remark 8. Conditions (22c) imply the continuity of Jh across interdomain bound¬ 
aries. 


Remark 9. The hybrid variable A h is the Lagrange multiplier associated with the 
continuity constraint of Jh • n and it can be regarded as an approximation of u at the 
grid nodes: 

Ah , -u\m h - 


7.3. Static condensation. Consider a generic element K e Th- From (22a), we 
can compute Jh\x by inverting the local flux mass matrix as 


(23a) 


J K ~ 


-1 

K 


[B'k + C K )u K + L k \ k 


1)2x1 


Substituting J K into (22b), we get Uh\K 
(23b) 


u K — H k g K + B k A k L k \ 


k 


1)1x1 


where 


H k := 


—B k Ak(Bk + C K ) + E_ 


K 


D lX 1 


Now, using (23b) in (23a), we can express J k in terms of the sole hybrid variable 
A k and of gx as 


(23c) 


J k — M k X k + b K 


where 

and 


Mk : = 


A-k (Bk + C k )Hk B K Ajf L k + A K l L 


k 


d2x2 


h K ■■ = 


Ak{Bk + C K )Hffg K 


d2x1 


The above procedure is well known as static condensation and corresponds to Gaussian 
elimination of the internal variables Jh\k and Uh\x in favor of A h\di< at the level of 
the element K G Th (see 0)- 
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7.4. The algebraic system. Enforcing the continuity of J\ at each internal node 
through (22c) yields the following N e i — 1 conditions 

(24a) Aj) = Jf‘( X l+1 ) i = 1, ..., M - 1 

which can be written in the form of a linear tridiagonal algebraic system for A, i.e., 


for i — 1,. 

.., M — 1: 


(24ba) 

\m£ Ai_! 

+ (M# - 

(24bb) 

[Ai = Uq, 

A n = W. 




Equations (24b) constitute a linear system of the form 
(24c) MX = f 

where the entries of the stiffness matrix M G are defined as: 

(24da) iA 

(24db) 


(24dc) 


Mi,i = M 2 f - M u 


Mij +1 — — M\ 


12 
M—l 


and the entries of the load vector f G 
expressions of these entries will be specified in Sect. |8.1 


are fi = bf l+1 — b% . The explicit 


8. The Stabilized DMH Method and the Discrete Maximum Principle 


If the solution of a given boundary value problem satisfies a CMP, then a properly 
designed approximation should behave in the same way. A numerical scheme that 
does not generate spurious global extrema in the interior of the computational 
domain is said to satisfy a discrete maximum principle (DMP). We assume that the 
solution u of the BVP (18) satisfies the a priori estimate ( |19h ) and we characterize the 
conditions under which the same estimate is satisfied also by the DMH approximation 
A h computed by solving the linear algebraic system (24c) associated with the problem. 

The following definitions turn out to be useful. 


Definition 2 (Inverse monotone matrix). An invertible square matrix A of size n is 
said to be inverse monotone if 

(25a) A -1 > 0 

the inequality being understood in the element-wise sense. 


A special class of monotone matrices is that introduced below. 

Definition 3 (M-matrix). An invertible square matrix A is an M-matrix if: 

• Aij <0 fori ^ j; 

• A is inverse-monotone. 


Given a function rjh G A/j iP , for any given function p : <9f 1 —> M, we denote by p* h 
the piecewise linear continuous interpolate of rjh over 7/, . 

Theorem 3 (Sufficient condition for DMP). Assume that matrix M is an M-matrix. 
Then \* h satisfies the DMP, i.e. 

(25b) 0 < A * h (x) < M, ^ ViGh. 












THERMO-FLUID-ELECTROCHEMICALLY DRIVEN ION FLOW 


19 


The following (necessary and sufficient) condition is useful to verify the property 
of being an M-matrix. 

Theorem 4 (Discrete comparison principle). Let A be an invertible matrix with 
non-positive off diagonal entries (Aij < 0 fori j). Then, A is an M-matrix if 
and only if there exists a positive vector e such that Ae >0 (in the component-wise 
sense), with at least one row index i* such that (He);* > 0. 


.1. The stabilized DMH method. The solution of the model BVP (18) may 


exhibit sharp boundary and/or internal layers according to the relative weight of 
the advective and reactive coefficients with respect to the diffusion term (see [[50]). 
Correspondingly, this is well known to reflect into numerical instability (spurious 
unphysical oscillations) if the local mesh size h is not sufficiently small. However, 
a smaller value of h means a larger size of the linear algebraic system and thus 
an increased computational effort. To avoid this inconvenience, we introduce in 


the DMH formulation (22) two specific approaches that ensure numerical stability 


without necessarily resorting to a small value of h : 

(51) : diagonal lumping of the local mass flux matrix Ar (for dominant reac¬ 
tion); 

(52) : addition of an artificial diffusion (for dominant convection). 

In particular, in the remainder of the section we show that, when (SI) and (S2) are 
used in conjunction, the stiffness matrix M. is an M-matrix irrespective of the value 


of h. This property guarantees that system (24c) is uniquely solvable and that the 


piecewise linear interpolate of the solution, X* h , satisfies the DMP and, in particular, 
the a-priori estimate (25b). We refer to m for a thorough discussion of continuous 
and discrete maximum principles enjoyed by the solution of singularly perturbed 
BVPs and by their corresponding finite element and finite difference approximations. 
For ease of presentation, we assume henceforth that problem coefficients D,v,c and 
g are constant positive quantities and that the mesh size is uniform and given by 
h = l/N el . 

8.1.1. (SI): Lumping of the local flux mass matrix. The stabilization approach to 
deal with the case where the reaction coefficient c dominates over the diffusion 


coefficient D in (18) consists of replacing the local flux mass matrix Ak in (23a) 


with a diagonal matrix A K such that: 


= J2(A K )ij = 

o =i 


(26a) 

(A k 

(26b) 

(A k 

The approach (26) 


i — 1,2 

i ^j- 


The approach (26) is referred to as mass lumping and is equivalent to using the 


trapezoidal quadrature rule to compute the integral 


where 


Using (26) in (22a), when q h = we get 
(26c) 


(■ A K )ij= 'ipjifidx i,j = 1,2 

J K 

P 1 (/l)= span {^, 1 / 2 } V/feT, 

Ur ~ Ai 


Ji = vu K - D h 


h/2 ’ 
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while, when rp, — f 2 , we similarly get 


(26d) 


J 2 = vu K - D, 


X 2 — uk 
1 h/2 ' 


Equations (26c) and (26c) express each local degree of freedom of the flux Jh\i< as 


the sum of an average advective flux and a finite difference approximation of the Fick 
diffusive flux flowing between the boundary of K and its center. The two relations 
are the DMH analogues of those of the dual-mixed stabilized method of 


8 .1.2. (S2): Artificial diffusion. The stabilization approach to deal with the case 
where the advection coefficient v dominates over the diffusion coefficient D in (18) 


consists of replacing the diffusion coefficient D with the following modified expression 
(cf. 0) 


(27a) 


where Ve ioc : = 


\v\h 

TD 


Dh D (1 + &(fPei oc )) 


is the local Peclet number and $ : M + —> M + is a suitable 


stabilization function to be chosen in such a way that: 
(27b) $(t) > 0 


(27c) 


lim $(t) = 0. 

t—► o+ 


The approach ( |27a ) is referred to as artificial diffusion. Two special choices of <f> are: 
• Upwind (UP) stabilization function 

$ UP (t) := t 

Scharfetter-Gummel (SG) stabilization function 

t SG (f) :=«-! + But). Bit) ■- * 


(27d) 


(27e) 


e — 1 


In the case where (27e) is adopted, the resulting stabilized DMH coincides 


with the classic SG exponentially fitted difference scheme (cf. 

The accuracy of the two stabilized DMH methods, as h —> 0, is O(h) if $ = <& UP 
while is 0(h 2 ) if <F = <h SG '. Thus, the SG formulation is far superior than the upwind 
method in the limit of a small grid size. However, for a large value of the Peclet 
number the two stabilized methods are practically equivalent in terms of the amount 
of added artificial diffusion (see also [33] and [25]). 


8.1.3. Matrix form of the stabilized DMH method. Using superscripts — and + to 
refer to quantities on the interval K~ = [xi-\,xf\ and Kf = [xi,Xi + 1 ], respectively, 
let us enforce continuity of J h at every interior node Xi [i — 2,..., N — 1) 

(28a) 


— T+ 

J 2 — d 1 . 


Then, combining equations (26c) and (26d) with condition (28a), we get 


(28b) - Dff + (l '~ + + 2A -°") 

v ; h 4D~ + c~h 2 


.2 D+ (v + ~^ r )(g+h 2 + 2\iD++ 2X i+1 D+) 

/ W?. I 


h 


4 D + + c + h 2 
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By inspection on (28b) the entries of the stiffness matrix A4 and of the load vector f 
become: 


(28ca) 


(28cb) 


M 


i,i —1 


(LVf 5 ) 


1 


1 + 


ch 2 
AD 


M 


l,l 


2D(1 + 4>) 

hT 


/ 

i + 


ch 2 \ 

2D 

ch 2 

AD/ 


(28cc) 


(28cd) 



fi 


gh 


1 + 


ch 2 ' 

AD 


D{ 1 + <I>)\ 
h ) 


1 


1 + 


ch 2 

AD 


Proposition 3 (Discrete Maximum Principle in the non stabilized case). Set $ = 0 
and assume that 


(28d) Vei oc < 1. 

Then, the stiffness matrix A4 is an M-matrix. 


Proof. Recalling that v > 0, it turns out that JAij -1 is < 0. Moreover, we have 
Mi : i > 0 and f > 0. By inspection on we see that if (28d) is satisfied then 

M. is an M-matrix and the DMP holds for the solution of (24c). □ 


Condition (28d) may be too restrictive in the choice of the mesh size so that the 


following result may be a convenient remedy. 


Proposition 4 (Discrete Maximum Principle in the stabilized case). When $ = <A> UP 
or <3> = < f> SG 7 the stabilized DMH method satisfies the DMP irrespective of the value 

of 'Peioc- 


Proof. It is easy to verify that 

E Mij > 0 

j=l,...,N 

for every interior node (i — 2,..., N — 1), and 

+ A4i j2 > 0 

+ M-n,n > 0 

at boundary nodes. Then, the conclusion immediately follows by applying Theorem |4} 

□ 


9. Experimental Validation of the DMH method 

In this section we numerically verify the convergence and stability properties of 
the DMH method. 
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N el 

| \U - U h | L 2 

||n 0 w — Uh\\L 2 

\\ u - K\\l 2 

|| M Xh ||oo,/l 

10 

4.77887e-01 

2.21777e-01 

1.30666e-01 

8.64963e-02 

20 

1.91869e-01 

6.25684e-02 

3.38644e-02 

2.22424e-02 

40 

8.63963e-02 

1.61248e-02 

8.54321e-03 

5.60052e-03 

80 

4.18183e-02 

4.06198e-03 

2.14066e-03 

1.40542e-03 

160 

2.07297e-02 

1.01743e-03 

5.35468e-04 

3.51510e-04 

320 

1.03422e-02 

2.54479e-04 

1.33886e-04 

8.78873e-05 

640 

5.16826e-03 

6.36272e-05 

3.34727e-05 

2.19724e-05 

1280 

2.58377e-03 

1.59073e-05 

8.36820e-06 

5.49311e-06 

2560 

1.29184e-03 

3.97683e-06 

2.09202e-06 

1.37325e-06 


TABLE 1. Various error norms for numerical solutions Uh and X* h as a 
function of the number of intervals N ei . No stabilization is adopted. 


N el 

a, 

1 

t-i 

to 

\\J ~ JhWH 1 

10 

3.28036e-01 

2.60239e+00 

20 

8.75295e-02 

1.30056e+00 

40 

2.16140e-02 

6.47536e-01 

80 

5.29217e-03 

3.23222e-01 

160 

1.30365e-03 

1.61529e-01 

320 

3.23124e-04 

8.07533e-02 

640 

8.04091e-05 

4.03752e-02 

1280 

2.00543e-05 

2.01874e-02 

2560 

5.00749e-06 

1.00937e-02 


TABLE 2. Error norms H 1 and L 2 for numerical solution Jh as a 
function of the number of intervals N e i. No stabilization is adopted. 


9.1. Convergence Analysis. We consider the model BVP (18) with with ho¬ 
mogeneous BCs, u(0) = u(L ) = 0, d, — 5, coefficients D — v — c — 1, and 
g(x ) = e~ x [:x 2 — (6 + L)x + 3L — 2] in such a way that the exact solution is the pair 

u(x) = xe~ x (L — x), J(x) = — e~ x 2x 2 + 2 x(L — 1) + L 


In the numerical approximation of the problem, the mesh size is uniform and equal 
to h = d/N e i, with N e i = [10, 20,40, 80,160,320, 640,1280, 2560] 7 . Table [l] illustrates 
the convergence history of the DMH method by reporting various error norms 
for u — Uh and u — X* h . The discrete maximum norm || • ||oo,a associated with the 
triangulation Th is defined for each continuous function rj(x) : [0, L\ aR as 


IWaOlloo,*. : = max |p(^)l- 

xieJV h 

Table [2] reports the error J — Jh in the L 2 -norm and in the II,\\v = H 1 - norm. 

The analysis of the asymptotic convergence orders that are predicted by Tables [T| 
and [2] show that the computed numerical solutions verify the following error estimates 

• ||w — Uh\\L 2 < Ch 

• |jn 0 M - Rh||L 2 < Ch 2 

• jju — X* h \\L 2 < Ch 2 

• |ju - A^ljo^h < Ch 2 
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• || J — J h \\ H i Ch 

• ||J- j h \\v <Ch 2 

The above results are in excellent agreement with the theoretical convergence rates 
predicted in the elliptic case in pjj and [5j [29j. 


9.2. Reaction-Dominated and Advective-Dominated Regimes. In this sec¬ 


tion we demonstrate the efficacy of the stabilization techniques proposed in Sect. 8.1 


in the study of two model problems, special instances of the BVP (18). For ease of 
presentation, we set L — 1 and we subdivide the computational domain into N e i = 10 
uniform intervals of size h = 1/N e i = 0.1. We also assume that the coefficients D : v : c 
and q are constant, with q = 1 and we take homogeneous boundary conditions, 
u(0) = u(l) = 0 (u = 0). 


9.2.1. Diffusion-reaction BVP: mass-lumping. In this case we have v = 0. A barrier 
function for u is <f(x) = 1/c. If the diffusion coefficient is small compared to the 
reaction term, e.g D = 10 -3 , c = 1, spurious (unphysical) oscillations arise near 


the boundaries, see Figure 2(a) In order to eliminate such oscillations, we adopt 


the mass-lumping stabilization procedure. The numerical solution A^ obtained with 


mass-lumping is illustrated in Figure 2(b) 




(a) No lumping 


(b) Lumping 


Figure 2. Exact and numerical solution of the diffusion-reaction 
BVP. Exact (green, solid line) and numerical solutions A h (blue line 
with bullets). 


9.2.2. Diffusion-advection BVP: artificial diffusion. In this case we have c = 0. A 
barrier function for u is f>(x) = x/v. If the diffusion coefficient is small compared to 
the advective term, e.g D = 5 ■ 10 -3 , v = 1, spurious (unphysical) oscillations arise 
in the neighborhood of x — 1 and propagate throughout the entire domain polluting 
the overall quality of the computed solution, see Figure |3(a) Numerical solutions A/,, 
obtained with UP and SG stabilization functions are illustrated in Figures |3(b) and 


3(c), respectively. In the case of the SG stabilization function, the computed A h is 


nodally exact. 

























24 


RICCARDO SACCO 1 , FABIO MANGANINI 2 , AND JOSEPH W. JEROME 3 



FIGURE 3. Exact and numerical solution of the diffusion-advection 
BVP. Exact (green, solid line) and numerical solutions A h (blue line 
with bullets). 


10. Simulation of biological channels 


In this concluding section we carry out a thorough validation of the vET and 
vTHD models in the study of two different biological channels: 1) the Gramicidin-A 
channel considered in £0J; and 2) the bipolar nanofluidic diode considered in [30]. 
The main focus of the simulations is on the current-voltage (IV) characteristics of 
the channel and on how the IV curves are affected by the boundary conditions, 
especially the temperature of the two bulk regions, T(0) = Tr and T(d ) = Tr. 
Also, we aim to investigate the electroosmosis effect, which is accounted for by the 
electrolyte fluid velocity v e . In all the reported computations we set N e i = 200 and 


use the SG stabilization function (27e), with the exception of a specific case discussed 


in Section |10.1| Moreover, in the case of the vTHD model we set the saturation 
velocity v sat equal to lOnrs” 1 . We also assume that the channel has a constant 
cross-sectional area equal to A , whose value may vary from channel to channel. For 
further information on the simulated channels, we suggest to consult [21] and the 
references cited therein. 


10.1. Gramicidin-A channel (ballistic diode). This channel is thoroughly an¬ 
alyzed in the work ra that is here used as a benchmark for the biophysical and 
numerical assessment of models and methods proposed in the present article. To allow 
comparison between the results of our models and those of ffl we set v e = 0. Unless 
otherwise specified, the subsequent figures refer to computations with the THD 
model, meaning that the outcome of the ET model would appear indistinguishable. 
The channel length d is equal to 2.5 nm, p>R is always set equal to 0 V and <pl = V app) 
V app being the externally applied voltage. 

10.1.1. Electrochemical variables. The permanent charge profile of the channel is 
illustrated in Fig. [I] (left). Because of the negative fixed charge the channel is highly 
selective to ion flow and attracts positive Na + ions while Cl~ ions are mainly repelled. 
This behavior is confirmed by the computed ion concentration profiles shown in 
Fig. [4] (right) in the case where a voltage drop (fiL — Pr = Va PP = 0.1 V is applied 
across the channel. We first observe that the cation concentration computed by our 
algorithm is in excellent agreement with that of m- We also notice the presence 
of strong internal layers in the Na + distribution which are effectively captured by 
the stabilized DMH discretization without introducing spurious oscillations. In 
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the mentioned reference the hydrodynamic equations are solved using a rather 
different numerical strategy based on the essentially non oscillatory finite differences 
with shock capturing introduced in [32, I3HJ- The role of the SG stabilization in 
the quality of the computed ion concentrations is well documented by the results 
shown in Fig. [5| In this graph, we see on the left side the cation and anion densities 
obtained by running the DMff discretization scheme without stabilization (<f> = 0) on 
an uniform grid of N e i = 19 elements. The maximum local Peclet number is in this 
case equal to 68.6288 so that the solution exhibits a markedly oscillatory behavior 
and fails to satisfy the DMP. Conversely, the adoption of the SG stabilization with 
the same number of elements produces the ion concentrations plotted on the right 
side of Fig. [5] The computed solution is strictly positive and the internal layers at 
the interfaces between the channel and highly doped regions are captured with the 
resolution allowed by the roughness of the grid size. 

Permanent charge and concentrations 


0 

S - 20 

-40 

5 

o 

'I -60 

a 

- 

S 

£ -80 

-100 

0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 

x(m) xiO- 9 x(m) x 10" 9 

FIGURE 4. Gramicidin-A channel. Left: permanent charge profile. 

Right: ion concentrations. 

The permanent charge profile and the consequent distribution of ions allows us 
to regard the Gramicidin-A channel as the biophysical analogue of an electronic 
"ballistic" diode of type p + -p-p + . This analogy is useful in the interpretation of 
simulation results and supports the idea that "channels are transistors alive" (cf. [T2j)- 
In this respect, the region of the channel x G [0, 0.5] nm corresponds to the Source 
terminal which has the role of emitting cations into the channel, while the region of 
the channel x G [2, 2.5] nm corresponds to the Drain terminal which has the role of 
collecting the cations that have travelled throughout the channel under the action of 
the thermo-electro-chemical forces. Electric potential and electric held profiles are 
shown in Fig. [6] (left) and Fig. [6] (right), respectively. Again, we observe that the 
electric potential and held computed by our algorithm are in excellent agreement 
with those of HD]. By inspection of the electric potential distribution, we see that the 
positive pole of the bio-electrical system is located at x = 0, nm while the negative 
pole is at x = 2.5 nm. Thus, the cations move from left to right, giving rise to 
accumulation of (hxed) negative charge at the entrance of the channel and positive 
(mobile) charge at the outlet of the channel. These accumulations correspond to the 
two strong peaks visible in the electric held distribution. 

To investigate the effect of different boundary values for the temperature at both 
entrance and outlet of the channel we let T(0) = T(d) range between 270 K and 
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Concentration profile without and with numerical stabilization 




Figure 5. Gramicidin-A channel. Concentration profile with a grid 
of 19 elements. When no numerical stabilization is adopted, spurious 
oscillations arise (left). The SG stabilization (right) leads to a stable 
solution. 


Potential and electric field 



FIGURE 6. Gramicidin-A channel. Left: electric potential. Right: 
electric field. 


330 K. Fig. [7] shows the IV curves for the THD model (left) and the ET model 
(right), respectively. The externally applied voltage V app ranges from 0 to 0.1 V. 
Temperature influence is visible in both models: the higher the bath temperature the 
lower the current. This agrees with the biophysical insight that an increase of bath 
temperature corresponds to a higher thermal energy for the ions and, consequently, 
to a higher energy dissipation because of frictional effects between particles and fluid. 
We notice also that the ET model predicts a slightly higher current than the THD 
model because this latter model accounts for convective energy dissipation within 
the fluid. 

10.1.2. Thermal variables. Cation temperature profiles when T(0) and T(d) are 
increased separately, are shown in Fig. [8j 
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IV curves when T( 0) = T(d ) 7 ^ 300 K 



Figure 7. Gramicidin-A channel. IV curves when T(0) = T(d) 7 ^ 
300 A' (left). THD model (left). ET model (right). 


Ion temperature when T( 0) or T(d) > 300 K 



FIGURE 8 . Gramicidin-A channel. Cation temperature profiles. 
T(0) > 300 K (left). T(d) > 300 K (right). 


Looking at the two families of distributions, we can distinguish five distinct 
subregions in each curve: the Source and Drain regions (reservoirs), the channel 
region and the two junctions separating Source and Drain from the ion channel. I 11 
the two reservoirs, the temperature profile is linear, because the cation concentration 
is almost constant and the electric field is very small. Then, ion temperature increases 
in the channel region according with the fact that particles are accelerated by the 
electric held from left to right, approximately for x > 0.8 11111 (cf. Fig. [ 6 j right). 
However, the effect of the electric held at the two junctions is quite different in the 
two sets of thermal boundary conditions. If T(0) > T(d), electric and thermal helds 
act in the same direction (from left to right) so that as T( 0 ) increases, the thermal 
how at the two channel boundaries correspondingly increases. If T(d) > T(0) electric 
and thermal forces act in opposite directions and this contributes to diminishing the 
thermal how at the two channel boundaries. It is interesting to notice that ion heating 
(i.e., the maximum value of ion temperature inside the channel) is considerably larger 
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in the case where T(d) > T(0), because in this condition ion acceleration due to 
the electric held greatly dominates over the thermal held along the channel (moving 
from left to right) so that the ion total energy increases and temperature gets larger. 
Then, once the ions are injected across the junction between channel and Drain 
they immediately thermalize to the local energy of the reservoir distribution and 
approximately cool down to the temperature enforced at x — d. 

Electrolyte fluid temperature when T(0) or T(d ) > 300 K 



FIGURE 9. Gramicidin-A channel. Electrolyte fluid temperature pro- 
flies. T(0) > 300/1 (left). T(d) > 300/1 (right). 


Temperature profiles of the electrolytic fluid are almost linear for both choices 
of the applied thermal drop, see Fig. [9j This behaviour is to be ascribed to the 
high value of thermal conductivity of the electrolyte fluid compared to that of the 
ion fluid. This makes the fluid behave as a perfect sink so that its heating is only 
passively driven by an external temperature gradient according to Fourier’s law (9f) 
(v e = 0 in this case). 


10.1.3. Ion velocity. From ion velocity profiles we can inspect how an externally 
applied temperature gradient may affect ion motion in the channel. 


Ion velocity when T(0) or T(d) > 300 K. THD model 




Figure 10. Gramicidin-A channel. Cation velocity computed by the 
THD model. T(0) > 300/1 (left), T(d ) > 300 K (right). 
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Ion velocity when T(0) or T(d) > 300 K. ET model 



Figure 11 . Gramicidin-A channel. Cation velocity computed by the 
ET model. T(0) > 300/1 (left), T(d ) > 300 K (right). 


In Fig. 10 (left) we see ions moving slower when T(0) > T(d) = 300 K. The 


opposite happens in Fig. 10 (right) where T(d ) > T(0) = 300 K. However, in this 


latter case, velocity variations are much smaller. Again, the different behavior of 
the ion velocity in the two sets of thermal conditions is related to the interplay 
between electric and thermal forces already discussed in Sect. 10.1.2. If T(0) > T(d) 


thermal diffusion is larger than in the case where T(d) > T(0), and increases as T(0) 
increases. This explains the maximum value of v n at T(0) = 300 K. If T(d) > T(0) 
ion motion is mainly driven by the electric held, which explains why the various 
curves are substantially insensitive to the increase of T(d). The ion velocity in the 
ET model exhibits a trend similar to that of the THD model, see Fig. [IT] even if the 
distribution of the peak value shows a much larger spread than in the case of the 
velocity computed by the THD model. This is to be ascribed to the larger thermal 
conductivity of the (thermally) unified system composed by anions, cations and fluid. 


IV curve when T(0) or T(d) > 300 K. THD model 



Figure 12. Gramicidin- A channel. I-V curves computed by the THD 
model. T(0) > 300 A' (left), T(d) > 300/1 (right). 
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IV curve when T(0) or T{d ) > 300 K. ET model 




Figure 13. Gramicidin- A channel. I-V curves computed by the ET 
model. T(0) > 300 A' (left), T(d) > 300 A' (right). 


The IV curves computed by the two models in the two distinct thermal sets of 
boundary conditions reflect what reported about ion velocities: a smaller ion velocity 
means a smaller current flowing in the channel. The externally applied voltage V app 
ranges from 0 to 0.1 V. Fig. 12 (left) shows the IV curves for the THD model when 
T(0) > T(d) = 300 K. Notice that the current assumes negative values when V app is 
close to zero and the temperature gradient due to thermal BCs is strong enough to 
counterbalance the affect of the applied voltage. Consistently with the ion velocity 
profiles of Fig. [lO] (left), the current flowing in the channel is lower than in the case 
of homogeneous thermal BCs. As T(0) increases, the I-V profile is shifted down 
with respect to the case of homogeneous BCs (T(0) = 300K solid black line). The 
case when T(d) > T(0) = 300K is shown in Fig. 12 (right). In this situation the 
IV relationships are shifted up, meaning that as T(d) increases, a higher current 
flows in the channel, just as predicted from the ion velocity profiles of Fig. 10 (right). 


Similar results and observations hold for the IV curves predicted by the ET model, 


see Figs. 13 


10.2. IV curves with velocity extended THD and ET models. In the sim¬ 
ulations conducted so far, we set v e = 0, which corresponds to neglecting the 
electroosmosis effect. We have also conducted simulations with v e varying in the 
range 0 -i-0.01 ms" 1 , and we found that it has no appreciable influence on the current 
flowing in the channel, so that no results are shown in this case. However, different 
conclusions are drawn in the channel presented in the following section. 


10.3. Bipolar nanofluidic diode. The nanofluidic channel investigated in the 
present section (whose detailed representation can be found in [TO], 2TJ) is called 
Bipolar (BP) nanofluidic diode because of its similarity with a p-n electronic diode: 
as a matter of fact, the BP channel can operate into two different states, depending 
on the sign of the applied voltage: (1) an ’open’ state (also called ’forward’ bias), 
characterized by high current flowing through the device, and (2) a ’closed’ state 
(also called ’reverse’ bias), with very little current flowing in the channel. 

10.3.1. Electrochemical variables. The input-output behavior of the BP channel is 
related to the surface permanent charge profile (cf. Fig. [I]). The channel has a 
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negative surface charge along half of its length, while the remaining part has positive 


surface charge (of the same magnitude), see Fig. 14 (left). Thus, anions and cations 
carry equal weight to channel behavior. This is confirmed by the symmetric spatial 
distribution of ion concentrations, both in the case of forward and reverse bias, see 
Fig. 14 (middle) and Fig. |l4| (right), respectively. Notice that ion concentration is 
much higher in forward bias than reverse bias. Electric potential and electric field 
profile are shown in Fig. [T5j One can see that, in the closed-state, carrier flow is 
inhibited by the potential barrier at the middle of the channel. Conversely, in the 
open-state, the potential drop enhances ion flow. The computed IV curves show the 
on-off trend just described. In particular, we see that if V app is negative the current 
is almost equal to zero, while it increases in a nonlinear manner if V app is positive (cf. 


Fig. 16). The predicted marked on-off function mode of the BP channel and the shape 


of the I-V curves are also in very good qualitative agreement with data reported in 
the study of heat biosensors in BHD ESI- This fact is an encouraging motivation 
to a further use and calibration of the computational models and methodologies 
proposed in the present article on biophysical novel applications. 

Permanent charge and concentrations 



FIGURE 14. BP channel. Permanent charge profile (left). Ion concen¬ 
trations: forward bias V app = +1V (middle), reverse bias V app = — IV 
(right). 


10.3.2. Temperature. The IV curves, when T(0) = T(d ) assumes lower or higher 


values than the reference value of 300 K, are shown in Fig. 16 for the THD model 
(left) and ET model (right), respectively. The externally applied voltage V app ranges 
from -1 to 1 V. The IV curves computed by the THD model follow the same trend, 
as a function of the bath temperature, as in Sect. 10. 1| (cf. Fig. [7]): the higher the 
bath temperature the lower the current flowing in the channel because of increased 
frictional effects, see Fig[l6] (left). The IV curves computed by the ET model in 
this type of channel differ remarkably from those computed by the THD model with 
respect to boundary temperature: the higher the bath temperature, the higher the 
current, see Fig 16 (right). We also observe that- the spread in the value of the 


maximum channel current as a function of bath temperature is much wider for the 
THD model than for the ET model. This trend was exactly the opposite in the 

as demonstrated by Fig. [7} We 


simulation of the Gramicidin-A channel of Sect. 10.1 


point out that the application to the BP channel of the theoretical prediction for 
channel current given by the ideal diode model (see, e.g., [2?]) would yield 

Q^app 


( 4 ) 


I = In 


exp 


KrT. 


- 1 


B -L sys / 
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Potential and electric field 




FIGURE 15. BP channel. Electric potential profile: forward (green) 
and reverse (red) bias. Electric field profiles: forward bias (green, left 
axis scale) and reverse (red, right axis scale). 


IV curves when T(0) = T(d ) 7 ^ 300 K 



Figure 16. BP diode. IV curves when T(0) = T(d) 7^ 300/1: THD 
model (left), ET model (right). 


where / is the total current flowing in the diode and / 0 is the saturation current. 
Thus, according to Q, an increase of channel system temperature turns out into a 
decrease in channel current, in accordance with the results computed by the THD 
model. Despite this preliminary indication in favor of the predictions of the THD 
model, we feel that the strongly different response of the two models when applied 
to different channel configurations certainly warrants further investigation and will 
be the object of a further step of our research activity. 

Also the temperature profiles of electrolytic fluid differ between the ET model and 
the THD model. The profiles from ET model have a non linear profile, while those 


from THD model are linear, see Fig. 17 (left) and Fig. 17 (right), respectively, when 
T(0) ranges from 270 to 330 K. 

10.3.3. Electroosmosis. In this type of channel the influence of the electrolytic fluid 


velocity is shown in Fig 18 (left) for the vTHD model and in Fig. 18 (right) for 
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T e (x ) profile when T(0) 7 ^ 300 K 



Figure 17. BP diode. Water temperature in the ET model (left) 
and THD model (right). 


the vET model. The selected range for the electrolyte fluid velocity (v e from 0 to 
0.01 m/s) is rather arbitrary since our model does not account, at the moment, for 
a self-consistent coupling between the vTHD system and the NS equations for the 
fluid. As the fluid velocity has positive sign, electrolytic particles flow from left to 
right. Accordingly, this additional translational force should enhance the movement 
of positive ions, while it should reduce that of negative ions, since the two ions move 
in opposite directions. From the IV curves shown, one can see that the additional 
driving force due to v e lowers the total current flowing in the channel in both models. 
Analogous results hold for negative values of v e (not shown). 

The temperature of the electrolytic fluid velocity for different values of v e is shown 
in Fig. [19} Temperature is only partially affected by v e . The increase in temperature 


predicted by the vET model (left) is about 3 K, which is quite remarkable, especially if 
compared to the results obtained with the vTHD model (middle), where temperature 
changes are much less significant, instead. The very different prediction of the vTHD 
formulation is related to the parameter v sat that strongly affects the value of the 
relaxation time parameter r Wv (cf. (lib)) in the energy balance equation (9c): the 
larger v sa t, the smaller t Wv . which corresponds to almost instantaneous restoration 
of equilibrium conditions in the electrolyte. Indeed, taking a smaller value of the 
saturation velocity, for example v sat = 1 m/s, leads to a result that is much closer to 
the thermal profile predicted by the vET model, see Fig. 19 (right). 


11. Conclusions and Research Perspectives 

In the present article we have proposed and numerically investigated a hierarchy 
of mathematical models for the simulation of thermal, fluid and electrochemical 
phenomena in biological transmembrane channels. The hierarchy is an extension of 
the classic Poisson-Nernst-Planck model for ion electrodiffusion and is conducted 
along the same lines of thought that have guided the development of the so-called 
hydrodynamic transport model in the analysis of semiconductor devices. To discretize 
the proposed models we have devised in the ID case a robust finite element dual-mixed 
hybridized method that ensures flux conservation, self-equilibrium and satisfaction of 
a positivity principle for ion concentrations and temperatures. The numerical scheme 
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IV curves when v e > Om/s 



FIGURE 18. BP diode. IV curve when v e ranges in 0 A- 0.01 m/s. 
vTHD model (left). vET model (right). 

T e (x) profile for different v e (vET vs. vTHD) 




Figure 19. BP diode. v e ranges in 0-E0.01ms _1 . vET model (left) 
differs greatly from vTHD at v sat = 10m/s (middle), while is quite 
close to vTHD at v sat = 1 m/s. 


has been thoroughly studied in several benchmark problems that demonstrate its 
accuracy and stability. An appropriate solution map is used to successively solve 
the nonlinear system of equations arising from model hierarchy and the resulting 
computational tool has been successfully calibrated and validated in the simulation 
of two realistic biological channels. Future developments of this study include: 

• time dependent simulations to describe the response of the channel to exter¬ 
nally applied stimuli; 

• self-consistent coupling of the hierarchy with the solution of the Navier-Stokes 
equations for the electrolyte fluid; 

• deeper investigation of the dependence of model predictions on biophysical 
parameters, for instance, saturation velocity that seems to play a critical role 
in determining the self-heating effect in the electrolyte fluid; 

• extension of the numerical scheme to 2D and 3D channel simulation. 
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